New code for equilibriums and quasiequilibrium initial data of compact objects 
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We present a new code, named COCAL - Compact Object CALculator, for the computation 
of equilibriums and quasiequilibrium initial data sets of single or binary compact objects of all 
kinds. In the COCAL code, those solutions are calculated on one or multiple spherical coordinate 
patches covering the initial hypersurface up to the asymptotic region. The numerical method used 
to solve field equations written in elliptic form is an adaptation of self-consistent field iterations in 
which Green's integral formula is computed using multipole expansions and standard finite difference 
schemes. We extended the method so that it can be used on a computational domain with excised 
regions for a black hole and a binary companion. Green's functions are constructed for various types 
of boundary conditions imposed at the surface of the excised regions for black holes. The numerical 
methods used in COCAL are chosen to make the code simpler than any other recent initial data codes, 
accepting the second order accuracy for the finite difference schemes. We perform convergence tests 
for time symmetric single black hole data on a single coordinate patch, and binary black hole data 
on multiple patches. Then, we apply the code to obtain spatially conformally flat binary black 
hole initial data using boundary conditions including the one based on the existence of equilibrium 
apparent horizons. 



I. INTRODUCTION 

In the last decades, simulation codes for compact ob- 
jects have been successfully developed in the field of 
numerical relativity, and various dynamical simulations 
have been performed includingbinary neutron stars and 
black holes inspirals to merger [l| , a massive core-collapse 
to a proto neutron star or a black hole (BH) formation 
, and black hole dynamics Q . Recent efforts on these 
subjects are, for example, to study more realistic situa- 
tions by including microphysics of the neutron stars (NS) 
[H, wider range of parameter space such as mass ratio 
and spins of binary black hole (BBH) mergers @ or BBH 
mergers in an ambient disk [7j. Accordingly, more re- 
alistic and accurate construction of initial data for such 
compact objects is required. 

Several researchers have developed methods for com- 
puting various types of initial data sets for those simula- 
tions Ic tilal. and equilibriums of rotating compact objects 
(see e.g. [16||). Many of such initial data codes are spe- 
cialized to a certain problem such as a single stationary 
and axisymmetric neutron star or BBH data on a confor- 
mally flat initial hypersurface. An exception is lorene 
[l7| , which is one of most used code for computing initial 
data for the merger simulations of binary neutron stars 
and black holes. The lorene code was originally devel- 
oped for computing rapidly rotating neutron stars, but 
has been extended to be capable of computing various 
kinds of equilibriums and quasiequilibrium initial data 
sets. 

In this paper, we introduce our project for develop- 
ing new codes for computing initial data of astrophysical 
compact objects, a single as well as binary compact ob- 
jects of all kinds, and present several tests for the new 
codes. Our aim is to develop a set of codes for comput- 
ing, on an initial hypersurface, a single neutron star (or a 



compact star such as a quark star), binary neutron stars 
and black holes, a central neutron star or black hole sur- 
rounded by a toroidal disk, and all these systems with 
magnetic fields. We call our new codes "COCAL" as the 
abbreviation for Compact Object CALculator 0. A note- 
worthy idea of the "cocal" project is to develop a code 
using less technical numerical methods than the recent 
initial data solvers with spectral methods [9 HlTi [l7| . Also 
the modules and subroutines of the fortran90 code are 
structured simply so that the code may be accessible by 
those who mastered introductory courses for programing. 
Such feature will help future developments to incorporate 
more complex physics in the code such as radiation, neu- 
trino radiation transfers, or realistic equations of state 
for the high density nuclear matter. 

The numerical method used in COCAL is based on 
Komatsu-Eriguchi-Hachisu (KEH) method for comput- 
ing equilibrium of a rotating neutron star [Lsl ]. In our 
previous works [HI, [l3| , we have extended KEH method 
for computing initial data for binary compact objects in 
quasiequilibriums. Among all, in the paper fl4| . we have 
introduced multiple spherical coordinate patches for com- 
puting binary compact objects. We improve the idea of 
the multiple patches in all aspects in the new COCAL 
code. In the paper [IH, we have presented convergence 
tests and solution sequences for rotating neutron star ini- 
tial data, which was calculated by the first version of 
COCAL. In this paper for introducing the COCAL code, 
we focus on the basic setup of the multiple spherical co- 
ordinate patches and the coordinate grids, the method 
of elliptic equation solver on the multiple patches, and 
convergence tests for binary black hole initial data. The 
paper is organized as follows: in SecUH we introduce an 



"Cocal" means "seagull" in Trieste dialect of Italian. 
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overview of the COCAL project, then coordinate setups, 
elliptic solver and other materials on numerical comput- 
ing. In Sec lIIII the results of convergence tests are pre- 
sented. In Sec lIVI solutions of BBH data on a conformally 
flat initial hypersurface are presented. We use geometric 
units with G = c = 1 throughout the paper. 



II. COCAL CODE 
A. Overview 

In the COCAL project, we aim to develop numerical 
codes for computing a single compact object as well as 
binary compact objects in (quasi)equilibrium using com- 
mon numerical method as much as possible. A plan for 
such codes also depends on how to formulate the prob- 
lem to solve such compact objects. Usually, a system 
of equations to describe equilibrium systems of compact 
objects involves a set of elliptic equations for the gravi- 
tational fields, and rclativistic hydrodynamical equations 
including Euler equation and the rest mass conservation 
equation, to which a stationary condition, either a time 
or a helical symmetry (l9j . is imposed. When the mag- 
netic field is present, elliptic equations for the electromag- 
netic fields arc added, and the equations for the fluid is 
replaced by magnctohydrodynamical (MHD) Euler equa- 
tion. Because the stationary Euler, or MHD-Euler, equa- 
tions are difficult to integrate numerically, a set of the 
first integrals in the form of algebraic equations, a suf- 
ficient condition for the stationary (MHD) Euler equa- 
tions, is derived and solved simultaneously with the field 
equations (see e.g. [20I l2lj). 

A choice of a numerical method is therefore made ac- 
cording to what kind of solver is used for solving the 
system of elliptic equations. The numerical method used 
in COCAL is based on KEH method for computing equi- 
libriums of rotating neutron stars [l8j . In this method, 
the elliptic equations are solved on spherical coordinates 
using the Green's formula iterativcly. This is done by 
separating the flat Laplacian or Hclmholtz operator on 
the variable to be solved for, then moving remaining (pos- 
sibly non-linear) terms to the source, and rewriting it in 
the integral form using Green's formula. Expanding the 
Green's function using spherical harmonics, the formula 
is integrated on the s phe rical coordinate grids numeri- 
cally (see, e.g. fl2l.fl3l.l22l]). The method is extended for 
computations of binary compact objects as discussed in 
this section. 

We choose simple finite difference formulas which are 
mostly second order accurate, and in some cases choose 
third or fourth order formulas only if they are necessary 
(see Sec. MI B 1[) . No symmetry, such as an equatorial 
plane symmetry, is assumed a priori on the 3D spherical 
computational domain. The COCAL code is written in 
Fortran 90 language, and runs with a few GB of memory 
for a model with a moderate resolution. 

We have developed basic subroutines for the COCAL 



code, including the coordinate grid setups for a single 
and multiple spherical coordinates, as well as the elliptic 
solvers for a single or binary compact objects, which we 
discuss in detail below. Mainly, two types of initial value 
formulation for Einstein's equation have been coded so 
far; one is assuming spatial conformal flatness (Isenberg- 
Wilson-Mathews formulation [23|, [13|), and the other 
non-conformal flatness (waveless formulation (l3l. l25j|). 

Also, the quadrupole formula to compute the gravita- 
tional wave amplitude and luminosity, and a Hclmholtz 
solver have been developed. We are in the phase to test 
all basic subroutines by computing simple test problems 
as well as known problems such as BBH initial data, or 
rotating neutron star solutions. In the next step, we will 
combine these developments and start computing new 
equilibriums and initial data sets such as helically sym- 
metric binary compact objects, or magnetized compact 
objects. 



B. Coordinate patches for binary systems 

We assume that the spacetime is M. is foliated by a 
family of spacelike hypersurface (£ t ) tg R, Ai = K x X 
parametrized by t 6 R. In the COCAL code, we solve 
fields on a initial hypersurface which may be sta- 
tionary (in equilibrium), or quasi-stationary (in quasi- 
equilibrium) . The initial hypersurface S t is covered by 
overlapping multiple spherical coordinate patches whose 
coordinates are denoted by (r, 9, (ft). Angular coordinates 
cover all directions (9, (ft) <E [0,n] x [0, 2tt] without any 
symmetry imposed. We also introduce Cartesian coor- 
dinates as a convenient reference frame in a standard 
manner, that is, to have the positive side of x-axis coin- 
cide with a (9, (ft) = (tt/2,0) line, that of y-axis with a 
(9, (ft) = (7r/2,7r/2) line, and that of z-axis with a 9 = 
line. 

In Fig. [TJ a schematic figure of three spherical coor- 
dinate patches whose coordinates are discretized in grid 
points is shown for the case of computing BH-NS binary 
systems by COCAL. Shown is the 2D section of the 3D hy- 
persurface, that may agree with the equatorial or merid- 
ional plane of the compact objects. Even though this 
may be the most complex setup for coordinate grids in 
COCAL, it is not technical at all compared to those of 
existing codes in which adaptive coordinates are used. 

For the computation of binary systems, two compact 
objects are placed at the centers of the two patches. 
We call these two patches the compact object coordi- 
nate patch (COCP), and the third patch the asymptotic 
region coordinate patch (ARCP). A domain of COCP is 
defined between two concentric spheres S a and Sb from 
which an interior of an another sphere S e is excised. 
Writing radii of S a , Sb, and S e as r a , rb, and r e re- 
spectively, we define spherical coordinates of COCP as 
(r, 9, (ft) £ [r a , x [0, tt] x [0, 2n], and locate the center of 
the excised sphere at (r, 9, (ft) = (d s , 7r/2, 0), that is on the 
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FIG. 1. A typical setup for multiple coordinate grid patches in the COCAL code for a BH-NS system. Left and right top patches 
are those for compact object coordinate patches (COCP) centered at each compact object (BH for the left, and NS for the 
right). The smallest circle with thick curve in COCP-BH is the sphere S a where the interior region is excised and certain 
BH boundary conditions are imposed. The ovals drawn in COCP-NS denote NS. Bottom patch is that for asymptotic region 
coordinate patch (ARCP), centered at the mass center of the system, and extends to asymptotics. The arrows represent maps 
of potentials between the multiple patches. Alternatively, the radius of COCP may be extended to asymptotics, instead of 
using ARCP. Note that the spheres S a , St, and S e of these coordinate patches are distinct ones on a spacelike hypersurface Et. 
The radius of each coordinate patch doesn't reflect the size used in actual computations. 
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positive side of the x-axis0- We introduce the excision of 
a domain interior of the sphere S e for computing binary 
systems, and elucidate its role in the following section. 
When a single and/or axisymmetric object is computed, 
the excision inside S e is not used. When BH is computed 
on COCP-BH with certain BH boundary conditions on 
S a (r a > 0), the region inside of S a is excised. When a NS 
or a puncture BH is calculated, the sphere S a of COCP is 
removed by setting r a = 0, so the radial coordinate cov- 
ers up to r = 0. A domain of ARCP is defined between 
two concentric spheres S a and Sb, and its spherical coor- 
dinates are defined as (r,6,(j>) e [r a ,r b ] x [0, n] x [0, 2ir]. 
When values of field potentials or other variables are com- 
municated from one patch to the other, those values on a 
certain sphere are mapped to a corresponding boundary 
sphere as indicated by arrows in Fig. [TJ 

Values of radii r a , r& , and r e for each of the coordinate 
patches used in actual computations will be summarized 
in Sec. Mil Typically, they are set as follows. For the 
case of using three patches as in Fig. [TJ the radius r a of 
the inner boundary S a of ARCP is taken large enough to 
be placed outside of the excised spheres S e for compact 
objects on COCP, but small compare to the size of the 
domain r b of the COCP. The outer boundary of ARCP, 
the radius r b of the sphere S b , is extended to the asymp- 
totic region when a field to be calculated behaves as a 
Coulomb type fall off, while it is truncated at the near 
zone when a radiation field is calculated. The center of 
ARCP is located at the center of mass of binary compact 
objects. Therefore typically, for compact objects with a 
mass M, r a = 0(M) (0 for NS), r b = O(100M), and 
r e = O(M) - O(10M) for COCP, and r a = O(10M), 
and r b = O(10 6 M) or larger for ARCP. 

As an another option for the choice of coordinate 
patch, the outer radius of each COCP r& may be ex- 
tended to asymptotics, say r b — 0(1O 6 M), and ARCP is 
removed. This option simplifies the code, but it can not 
be used when a radiation field is computed by solving 
Helmholtz equation. We present tests for the Helmholtz 
solver in a separate paper. 



C. Elliptic equation solver 

As mentioned earlier, the formulation for computing 
(quasi-)equilibrium configurations of compact objects re- 
sults in a coupled system of elliptic equations, cither Pois- 
son or Helmholtz equations with non-linear source terms, 
coupled with algebraic equations. The numerical method 
used in COCAL to solve such a system of equations is 
an extension of KEH method which is an application of 
self-consistent field method for computing equilibriums of 
self-gravitating fluids to general relativistic stars . A 



2 The positive side of x-axis of COCP-NS in Fig.[T]is pointing from 
the center of the coordinate grids toward left. 



distinctive feature of these methods is the use of Green's 
formula for an elliptic equation solver. We have intro- 
duced in previous papers our implementation of KEH 
method to compute binary neutron stars and black holes 

EH1. 

In the COCAL code, we have made a major change in 
the choice of coordinate patch, and accordingly in the el- 
liptic equation solver. Our new implementation is better 
in all aspects for computing binary compact objects than 
our previous ones. We come back to discuss this point 
after we introduce the elliptic equation solver in COCAL. 

In solving each field equation, we separate out a flat 
Laplacian or Helmholtz operator C, and write it with a 
non-linear source S, 



= S, 



(1) 



on an initial slice £j, where $ represents metric poten- 
tials. For the case of Laplacian L = A, using the Green's 
function without boundary G(x,x') = l/\x — x'\ that 
satisfies 



AG(x,x') = -4ttS(x-x'), 
Green's identity is obtained by 



(2) 



1 

47T 



G{x,x')S(x')d 3 x' 



-in 



+ — / [G(x, x')Y a $(x') - $(x')V a G(x, x')} dS' a . (3) 



dv 



where V is the domain of integration, x, x' £ V C S , 
and dV is its boundary. For the case of BH-NS system 
shown in Fig. [TJ the boundary of COCP-BH becomes 
dV = S a U S b U S e , that of COCP-NS dV = S b U S e , and 
that of ARCP dV = S a U S b . For the evaluation of the 
integrals in Eq.Q, a multipole expansion of G(x,x') in 
associated Legcndrc functions on the spherical coordinate 
is used; 



G(x,x') = 



1 



(£-m)\ 



\ x ~ x '\ tr ^ ^ +m y- 

xP £ m (cosd) P e m (cosd') cosm(^ - cp'), (4) 
where the radial Green's function ge(r,r') is defined by 



t(r,r') 



J+i ■ 



(5) 



with r> := sup{r, r'}, r< := inf {r, r'}, and the coeffi- 
cients e m are equal to eo = 1, and e m = 2 for m > 1. 

Eq.([3]) is an integral identity but is not a solution of 
Eq.CQ) in a sense that both of $ and its derivative n a V a $ 
can not be specified freely. Eq. ((3j> can be used to com- 
pute a potential over V, only if correct values of $ and 
n a V Q $ are known at the boundary dV. Here, n a is an 
outward normal to dV . Therefore, as it is commonly 
found in standard text books for clcctromagnctism [28| , a 
homogeneous function F(x, x') for the Laplacian is added 
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to evaluate the Green's function that satisfies the bound- 
ary condition at dV. For example, a Green's function 
G(x, x') + F(x, x 1 ) — at dV is used to impose Dirichlct 
boundary condition. In our previous paper [l4j , we have 
developed an elliptic equation solver on multiple coordi- 
nate patches that uses such Green's functions and solve 
Eq.© by iteration. 

A construction of such a Green's function that satisfies 
a boundary condition is, however, possible only when a 
certain specific geometry of the domain of computation 
is adapted to the coordinate systems. In the present case 
for COCP of the COCAL code in Fig. [TJ a Green's func- 
tion that satisfy boundary conditions at S a and Sb may 
not be derived in a practical form of equation, because 
we excised the region inside of the sphere S e . To im- 
pose boundary conditions at S a and Sb, we introduce a 
homogeneous solution xi x ), an d write a formal solution 
as 

$(ar) = x(x) + $int(», ( 6 ) 
where $int is equal to the right hand side of Eq.Q; 



$iNT(aO = -~r 

47T 



1 

47T 



G(x,x')S(x')d 3 x' 
[G(x, x')V' a 4>(x') - $(.t') V a G(x, x')] dS' a . (7) 



The homogeneous solution is computed so that the po- 
tential $ satisfies the boundary conditions, which are ei- 
ther one of Dirichlet, Neumann, or Robin boundary con- 
ditions at the boundary spheres S a or 5b; 



Dirichlet : $bc = /d 
Neumann : n a V a $BC = /n 



(8) 
(9) 



Robin 



f . 



n a V Q $BC + jV n°$ B C = /r (10) 



where /d, /n, and /r are given functions on the spheres 
So or Sb- Formulas for x( x ) are derived by using a Legen- 
dre expansion in an usual manner as shown in Appendix. 
Noticing x( x ) = &(x) — ^int^) the formulas for x( x ) 
can be written analogously to the surface integral terms 
of Green's formula, but with a different kernel function 
G BC 

x( x ) = ±- / [G BC (x, x')V a ($ BC - *mr)(x / ) 
47r Js a us b 

-(*BC - $INT)(x')V a G BC (x,x')] dS' a , (II) 

The function G BC (x, x') is expanded in terms of the as- 
sociated Legendre functions 



G BC {x,x') 



(£-m)\ 



71 {t + m)\ 

xPe m (cos6)P e m (cos6')cosm(<p-<p'), (12) 

where the radial function gf c (r, r') is chosen according 
to the type of boundary conditions used. We derive such 



G BC (x, 
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) Boundary S a 


Boundary Sb 


G NB (x 
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) None 


None 


G DD (x 


X 


) Dirichlet 


Dirichlet 


G m (x 


X 


) Neumann 


Dirichlet 


G DN (x 


X 


) Dirichlet 


Neumann 


G NN (x 


X 


) Neumann 


Neumann 


G RD (x 


X 


) Robin 


Dirichlet 


G DR (x 


X 


) Dirichlet 


Robin 



TABLE I. List of Green's function available in the COCAL 
code. The second and third columns correspond to the types 
of boundary conditions imposed on the boundary spheres S a 
and St, respectively. The case with no boundary condition is 
denoted by None. 



radial functions used in the corresponding surface inte- 
grals for various cases of boundary conditions as listed in 
Table |TJ Concrete forms of these functions are presented 
in Appendix [Bj 



D. Iteration procedure 

The final solution will be obtained from the iteration 
of Eq. ([5]), with Eqs. (J7J and (fTTj). where explicit form 
of Eq. (jlip for x( x ) depends on the boundary condition, 
for example, Eq. ([A"9]) or (|A12|) . 

We summarize the n th step of the Poisson solver in the 
COCAL code as follows: 

1) Compute the volume source term as well 
as the surface source terms on all possible surfaces 
S a , Sb, S e . 

2) Compute the volume integral and the surface inte- 
gral at S e for obtaining $int(^) from Eq. (|7J|. 

3) Compute the effective source for the integral on S a 
and Sb- For Dirichlet boundary condition it will be 
$bc - $int, while for Neumann - . 

4) Compute the surface integrals at S a and Sb for ob- 
taining x(x) according to Eq. ([lip using the appro- 
priate function G for the boundary conditions of 
the problem. 

5) Add the results from steps 2) and 4) to obtain $(x) 
from Eq. $6§. 

6) Update $^™^ according to 

$ {n \x) := S(x) + (1 - c)& n -V{x) 
where 0.1 < c < 0.4. 

7) Check if 

l$(n) _ $(«-!) I 
2^^ ', < e n 



|$(")| + |$("-1)| 
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for all points of the grids, where e c 



10" 



10" 

-7 



is taken in typical computations, and e c = 10 in 
this paper. If yes exit. If no go back to step 1). 

Here, intermediate variables during an iteration step are 
denoted with a hat as $int, Xj an d The above iter- 
ation procedure is applied to each coordinate patch one 
after other. In step 1), the sources of the surface terms 
are computed either from boundary conditions to be im- 
posed on the surface, or from data of corresponding sur- 
face on the other patch (see, Fig. [1] how the potentials 
are transferred from a boundary surface to the other). 
Several different iteration schemes are possible for solv- 
ing a set of elliptic equations for more than one variable. 
As long as we experimented, a convergence of the iter- 
ation does not depend on the order of computing those 
variables at each iteration step. 

We will see this elliptic equation solver produce accu- 
rate solutions for test problems of binary black hole data. 
Two comments on the elliptic equation solver are made 
here. Although, $int in Eq- © involves surface integrals 
on all S a , Sb, and S e , those on S a and Sb are not included 
in $int in an actual computation. Those computations 
are redundant because the homogeneous solution xi x ) 1S 
determined again from the surface integrals on S a and 
Sb as in Eq. ([TT]). So far, we do not plan to develop el- 
liptic solvers for vector (tensor) fields in which Green's 
functions are expanded in vector (tensor) spherical har- 
monics. Instead, we write the Cartesian components of 
vector or tensor equations, and solve each components as 
scalar equations on spherical grids for simplicity. We will 
see an example in Sec. |IV] (see also, [HI). 



E. Grid spacing 

We apply finite difference scheme to solve the system 
of equations for compact objects on the spherical domain 
introduced in Sec III Bl (see, Fig. [TJ. Spherical coordi- 
nates (r, 9, 4>) for COCP and ARCP are bounded by two 
concentric spheres S a and Sb of radius r a and r^, respec- 
tively, with the possible excision of a sphere S e of radius 
r e inside COCP. The origin of the radial coordinate r is 
placed at the common center of S a and Sb, where the 
compact object is placed for the case of COCP. Excised 
sphere for a binary companion S e is always positioned at 
a positive value on the x-axis at a distance d s from the 
origin. Clearly r a < d s — r e . For neutron star calcula- 
tions the sphere S a is absent and the coordinate system 
extends from r = to rb- 

In the COCAL code, the spacing of all coordinate grid 
points (n, 9j, 4>k) with i = 0, ■ ■ • , N r , j = 0, • ■ ■ ,N$, and 
k = 0, • • • , Nj,, are freely specifiable. However, in (6, <j>) 
directions, uniform grids are recommended to resolve the 
trigonometric and associated Legendre functions used in 
the elliptic equation solver, as well as a structure of com- 
pact objects evenly. That is, we set the grid interval in 



these directions as 

A8j = 6j - = A6 = 



A<p k = 



TT 



(13) 
(14) 



The grid spacing in the radial direction r is usually 
constructed on one hand to resolve the vicinity of the 
compact object with finer grid spacings, and on the other 
hand to extend to asypmtotics using increasingly sparse 
spacingf0. The setup for radial grids of COCP in the 
present computation is illustrated in Fig. [2j The grid 
is composed from three regions I, II, and III. For the 
case with r a < 1, the region I is set by r G [r a , 1] the 
region II by r £ [l,r c ], and region III r e [r c ,r*(,]. We 
introduce grid numbers N*, A r ™ which correspond to the 
numbers of intervals in regions I and I+II, respectively. 
We introduce a standard grid spacin g Ar as Ar = l/N s r . 
For the case with r a < 1, the grid intervals, Ar; := r, — 
ri-i, are defined by 

Ar l+l = hAr h for i = 1, • • • , N f r - 1 (15) 
Ar, = Ar, for i = A^ f , • ■ ■ , N? (16) 
Ar l+l = fcAri, for i = N™, ■ ■ ■ , N r - 1 (17) 

which correspond to regions I, II, and III in Fig. [51 respec- 
tively, where ratios h(< 1) and fc(> 1) are respectively 
determined from relations 



1 



n 



Ar- 



1 



= A: 



1-h ' 

' k- 1 ' 



(18) 
(19) 



For the case with r a > 1, which is mostly for ARCP, the 
grid intervals, Ar,, are defined by 



A^ = Ar for i = 1, • ■ ■ ,A^ r m , 
Ar i+ i = kAn for i = iV r m , • • ■ , N r 



(20) 
(21) 



where the ratio k is determined from Eq. (fl9| . Parame- 
ters for the grid setup are listed in Table QTJ 



F. Finite differences and multipole expansion 

Approximations made in our numerical method are a 
truncation of the series of Legendre expansion at a finite 
order of multipole, and an evaluation of a solution on 
discretized grids - the finite differencing. The accuracy 
of the code is, therefore, determined from finite difference 
formulas to be used, the number of grid points and their 
spacings, and the number of multipoles being included. 

In the COCAL code, we use second order mid-point 
rule for numerical integrations and differentiations, along 



3 For the case of solving Helmholtz equation the region extends to 
a near zone (a size of a several wavelengths of a dominant mode) . 
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FIG. 2. Radial coordinate grids for COCP for the case of regions for BH and binary companion being excised. The radial 
coordinate grids corresponds to those of lowest resolutions Al, Bl, Dl and Fl in Tables ITTT1 and ITVl 



r a ■ Radial coordinate where the radial grids start. 
rt : Radial coordinate where the radial grids end. 
r c : Radial coordinate between r a and rt where 

the radial grid spacing changes from 

equidistant to non-equidistant. 
r e : Radius of the excised sphere. 
N r : Number of intervals Ar^ in r G [r a ,r;,]. 
A'"' : Number of intervals Ar; in r G [r a , 1] for r a < 1. 

or r G [r a ,r a + 1] for r a > 1. 
: Number of intervals Ar; in r G [r a ,r c ]. 
Ng : Number of intervals A9j in 9 G [0, 7r] . 
N$ : Number of intervals Acf>k in G [0, 2ir]. 
d : Coordinate distance between the center of S a (r = 0) 

and the center of mass. 
d s : Coordinate distance between the center of S a (r = 0) 

and the center of S e ■ 
L : Order of included multipoles. 

TABLE II. Summary of grid parameters. The radius r e is 
defined only for COCP. 

with second order linear interpolation rule. In the elliptic 
equation solver ((6]) , the source terms are evaluated at the 
mid-points 

, a a \ f r i + r i-i Sj+Sj-i 0fe + 0fe-i\ 
(r i _ h ,B j _i,<f> k _i)= ^ , , j 



and integrated with the weights AriAOj A(f>k (other than 
a Jacobian). The mid-point rule has a few advantages. 
The second order accuracy of the mid-point rule for a 
quadrature formula is maintained even with a disconti- 
nuity of the derivative of Green's function for a volume in- 
tegral at r = r' . It may be possible to derive a higher or- 
der quadrature formula for numerically integrating such 
functions, but for instance, a Simpson rule doesn't guar- 
antee fourth order accuracy at grid points with i being 
odd integers. Also, an excision of a region inside of a 
sphere S e for a binary companion on COCP complicates 
a derivation of a higher order quadrature formula which 
maintain the degree of precision near the sphere. Because 
of the simplicity of the mid-point rule, it is not difficult 
to modify the weights for an integration to maintain the 
accuracy. Another advantage of the mid-point rule is to 
avoid the coordinate singularities of the spherical coordi- 
nates. 

In some cases, we also use third or higher order finite 
difference formula for the numerical differentiations. Es- 
pecially, it is found that it is necessary to use the third 
order finite difference formula for the radial derivatives 
to maintain second order convergence of the field near 
the BH (see, Sec lIII B Tj) . Interpolations of scalar func- 
tions from the grid points to the mid-points are done 
using a second order linear interpolation formula. We 
often need to interpolate a function from one coordinate 
patch to the other, such as to compute source term at 
S e of COCP. In such case, the functions are interpolated 
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Type 


r a 


Tb 


r c 


A r f A r m 


N r 


Ng 




L 


Al 


0.2 


10 4 


1.25 


16 


20 


48 


48 


96 


12 


A2 


0.2 


10 4 


1.25 


32 


40 


96 


48 


96 


12 


A3 


0.2 


10 4 


1.25 


64 


80 


192 


48 


96 


12 


A4 


0.2 


10 4 


1.25 


128 


160 


384 


48 


96 


12 


Bl 


0.2 


10 4 


1.25 


16 


20 


48 


24 


96 


12 


B2 


0.2 


10 4 


1.25 


16 


20 


48 


48 


96 


12 


B3 


0.2 


10 4 


1.25 


16 


20 


48 


96 


96 


12 


B4 


0.2 


10 4 


1.25 


16 


20 


48 


192 


96 


12 


Dl 


0.2 


10 4 


1.25 


16 


20 


48 


24 


24 


12 


D2 


0.2 


10 4 


1.25 


32 


40 


96 


48 


48 


12 


D3 


0.2 


10 4 


1.25 


64 


80 


192 


96 


96 


12 


D4 


0.2 


10 4 


1.25 


128 


160 


384 


192 


192 


12 



TABLE III. Grid parameters used in convergence tests for 
a single BH and equal mass BBH data solved on a single 
coordinate patch, COCP. 



using fourth order Lagrange formula. For example when 
the surface integral at the excised sphere is computed we 
need the potential and its derivative at point x' on S e as 
seen from the center of S e . Theses values are taken by 
interpolating the nearby 64 points of the other coordi- 
nate system. Some examples of finite difference formulas 
frequently used in the COCAL code is summarized in Ap- 
pendix [Cj 

As discussed in [l4[, the excised region is introduced 
to improve the resolution in angular directions, and ac- 
cordingly to reduce the number of multipoles to resolve a 
companion object. Without the excised region, the size 
of the companion object itself has to be resolved by an- 
gular grids, while in our setup, it is enough to resolve the 
size of the excised region, which is usually taken as large 
as the half of separation ~ d s /2. Then an angle to be 
resolved can be always about ~ 2arcsinl/2 = 7r/3. Note 
that although it is, in principle, possible to excise S e with 
a different radius from each COCP, and it is allowed in 
COCAL, it is more practical to have the same size of ex- 
cised region for the same reason above. To summarize, 
the angular resolution of a COCP is determined from 
the degree of accuracy to resolve the deformation of the 
compact objects centered at the patch, and to resolve the 
size of excised sphere, ~ n/3. The angular resolution of 
ARCP depends just on how many multipoles one wish to 
keep in the near zone to asymptotics. For both COCP 
and ARCP, the number of Legcndrc expansions are in 
the range £ ~ 10 — 16 for computing binary systems. 

Legendre P^" 1 may have t zero crossings in 9 6 [0, 7r), 
and sinm</> or cosm0 have 2m zeros in cj> £ [0, 27r). The 
number of grid points along the angular coordinates has 
to be enough to resolve these multipoles with maximum 
I or m, say 4 times more than the number of zeros. 



III. CODE TESTS 
A. A toy problem for black holes 

Convergence tests for a time symmetric BH and BBH 
data are performed to check the numerical method of CO- 
CAL presented in Sec|TTJ We assume the spacetime M is 
foliated by a family of spacelike hypcrsurfaces (E t ) tS R, 
A4 = M. x E parametrized by t G R. To obtain simple 
black hole solutions on E t , we assume time symmetric 
initial data, the extrinsic curvature K a b on E t vanishes, 
or in other words, assume the line element at the neigh- 
borhood of E( 



ds 2 



-a 2 dt 2 + ^p 4 f ij dx 1 dx j , 



(22) 



where is the flat spatial metric. Decomposing Ein- 
stein's equation G Q( g = with respect to the foliation 
using hypersurface normal n a to Et, and the projection 
tensor j ab = g al3 + n a rfi to it, we write the Hamilto- 
nian constraint G a ^n a n^ = 0, and a combination of the 
spatial trace of Einstein's equation and the constraint 
G a p(l a0 + \n a nP) = 0, as 



V 2 ^ = and V 2 (aV) = 



(23) 



These equations have solutions, which correspond to the 
Schwarzschild metric in isotropic coordinates for a single 
BH. For a two BH case, a BBH solution is given by Brill- 
Lindquist [29| 

Mi M 2 Mi M 2 , . 

1> = 1 + -± + and aij = 1 - -i - -1, 24 
2ri 2r 2 2r x 2r 2 

where subscripts 1 and 2 corresponds to those of the first 
and second BH; r\ and r 2 are distances from the first 
and second BH, respectively, and Mi and M 2 are mass 
parameters. The coordinates T\ and r 2 are written in 
terms of each other, for example in the first coordinate 
system of COCP, the radial coordinate are 



ri 



r and r 2 



yr 2 



d 2 — 2rd s sin 9 cos < 



where 9, </> the angular spherical coordinates of the first 
coordinate system, and 1 ■n- 2 for the second COCP. On 
the third coordinate system of ARCP, 



'2 



• 2 + d\ — 2rd\ sin 9 cos <j>, 
■ 2 + d\ — 2rd 2 sin 9 cos <f>, 



where d\ and d 2 are the distance between the center of 
ARCP and one of two COCP, and hence d s = d\ + d 2 . 

Instead of solving two Laplace equations Eq. (|23|) . we 
write a equation for a with a source on the whole domain 
of E t 



V> = and V 2 a 



(25) 
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In an actual computation, the BH centered at the COCP 
is excised at the radii r a of S a , and the binary companion 
is excised at the radii r e of S e which is centered at x = 
d s . Boundary conditions for these elliptic equations at 
S a are taken from analytic solutions (|24|) when Dirichlct 
boundary conditions are imposed. Neumann boundary 
conditions can be imposed with the use of 

dri 2r\ 2r| dr% 
da 1 (My M 2 dr 2 

For the outer boundary conditions, we choose Dirichlct 
boundary conditions whose data is taken from the ana- 
lytic solution Eq. (f2~4")l in all tests in this section. When 
the third parch, ARCP, is not used, Dirichlet data is im- 
posed on Sb (r — rfc) of COCP, while ARCP is used as 
in Fig. [TJ Dirichlet data is imposed only at Sb of ARCP. 

B. Convergence tests 

Convergence tests are performed to examine that the 
code produces solutions with an expected order of finite 
difference errors, and to find experimentally an (almost) 
optimally balanced set of resolutions for each coordinate 
grid (ri,9j,(j)k) which is not over resolved in one coordi- 
nate direction so as not to waste the computational re- 
sources. We find, from convergence tests for a single BH 
solution, it is necessary to use third order finite differ- 
ence formula for a radial derivative in the volume source 
terms in Eq. ([7]). We also find an optimally balanced 
resolutions between and 9j grids. From convergence 
tests for BBH data, we find appropriate resolution for 
4>i direction, and number of multipoles. Results for the 
convergence tests are discussed in this section. 

1. Single BH 

For the first test, we compute a single BH solution with 
mass parameter Mi = 2r a = 0.4 (and M 2 = 0). Eqs. (f2"5|) 
are solved on a single patch with a single excision region 
interior of S a for the BH. In Fig. [3l a fractional error of 
the lapse a 



8a 




*-* ^cxact 


a 




O^exact 



(28) 



is plotted along the x-axis for different resolutions in the 
radial coordinate grids r» (top and bottom left panels), in 
the zenith angle grids 9j (top right panel), and in all grids 
(bottom right panel). These resolutions are tabulated in 
Table |Hll and are indicated by A1-A4, B1-B4, and Dl- 
D4, respectively. In the set A1-A4, the radial resolution 
Ar is doubled, in B1-B4, the zenith angle resolution A9 is 
doubled, and in D1-D4, the resolutions in all directions 
are doubled at each level. Another difference in these 



results is the order of the finite difference formula used 
to compute a radial derivative in the volume source term 
in Eq. ([7]), where the second order (mid-point) formula 
is used in the top left and right panels, and third order 
(Lagrange) formula is used in the bottom left and right 
panels. It is noticed from the top left panel in FigJ3]that 
the error does not decrease as 0(Ar 2 ) when the number 
of radial grid points is increased as the parameter sets 
A1-A4 in Table ITO1 even for such a spherically symmetric 
solution. 

It appears that there are two reasons for that. In the 
top right panel, a convergence test is performed chang- 
ing the number of grid points in zenith angle 9j as the 
parameter sets B1-B4. It shows an improvement of the 
accuracy in 0(A9 2 ) in the larger radius r > 100, that 
is the error in this region is dominated by the finite dif- 
ferencing in 9 direction. However, the accuracy near the 
BH is not improved in both tests. In the bottom left 
panel of Fig. |3l the same convergence test as the top left 
panel is performed, but the finite difference formula for 
the radial derivatives is replaced by that of third order 
Lagrange formula 0(Ar 3 ). It shows that it is necessary 
to set the order of finite difference formula for the radial 
derivative as 0(Ar 3 ) to see 0(Ar 2 ) accuracy near the 
BH. This 0(Ar 2 ) error must be due to the mid-point 
rule used in the numerical integration. While, the error 
in the larger radius docs not decrease in the outer region 
with the radius r > 100. Finally, as shown in the bot- 
tom right panel of Fig J3] the error decreases in the second 
order for the set D1-D4, in which the third order finite 
difference formula is used for the radial derivatives. Con- 
vergence test are done also for increasing grid points in 
<j) directions <pk and also for changing the number of Lcg- 
endre expansion as L = 4 — 10, but they didn't change 
the results for such spherically symmetric BH test. 



2. Equal mass BBH computed with a single patch 

In FigfJJ results of convergence tests for equal mass 
BBH data arc plotted. In this test, we used only a sin- 
gle patch shown in FiglSl where the potential at the ra- 
dius r = r e rotated by tt in <f> coordinate is mapped 
to the excision sphere S e on the same patch to com- 
pute the equal mass data when the elliptic equations 
are solved. This amounts to impose the 7r-rotation sym- 
metry about the center of mass which is located at 
(r, 9, 4>) = (d s /2, 7r/2, 0). In this test, the number of grid 
points is chosen as D1-D4 in Table IIII1 separation be- 
tween the coordinate centers of two BH (a distance be- 
tween the centers of S a to S e ) is set as d s = 2.5, the 
excision radius of the BH r a = 0.2, the excision radius of 
the binary companion r e = 1.125, and the mass param- 
eters Mi = M 2 = 2r a . Hereafter, the third order finite 
difference formula is always used for computing the radial 
derivatives as discussed in Sec. IIIIB 11 

In the top panel of Fig.lU fractional errors of the lapse 
|<5a/a| defined in Eq. (|2"5|) are plotted along the x-axis 
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FIG. 3. Plots of fractional errors in the lapse 8a/a along the positive a;-axis. Plots show the errors with changing the number 
of radial grid points n as A1-A4 (top left panel), zenith grid points 6i as B1-B4 (top right panel), radial grid points n as A1-A4 
(bottom left panel), and all grid points (ri,6j,^>k) as D1-D4 (bottom right panel). In each panel, solid (red), long dashed 
(green), dashed (blue), and dotted (magenta) lines are in order from lower to higher resolutions. In the top panels, second 
order finite difference formula is used for calculating radial derivatives, while the third order formula is used in the bottom 
panels. 



near the BH. Because of the excision of the interior of 
the sphere S e and of the use of Legendre expansion in the 
elliptic solver, a certain modulation is seen in the errors. 
Therefore, hereafter we show fractional errors averaged 
over the number of (6j, (f>k) grids points at a radius r = ri 
defined by 



5 a 
a 



— r 



^cxact 



Q^cxact 



(29) 



where writing a grid point (ri,6j,<f>k) by p, we define a 
set Q t by Q t := {(r,, 0j, 4>k) \ V 6 V \ S? and n = const} 
where 5™ is an interior domain of S e . Then, is the 

number of points included in ft. 

The averaged fractional errors for the lapse are plotted 
along the radial coordinate r in the middle panel of Fig. [3] 
As expected, second order convergence is observed when 
the grid points arc increased as D1-D4. In the figure, it 
is seen that a couple of grid points at the vicinity of BH 
boundary S a have (averaged) errors as large as ~ 1% even 



for the highest resolution D4. This is due to our choice 
of the boundary r a = M/2 as same as the single BH test 
in the previous section. With this choice, the value of 
a becomes negative at the BH excision radius r = r a . 
Hence, the fractional error diverges at radii r (depend- 
ing on 9, and <p) where a crosses zero, even though the 
grid points are slightly off from the zeros. In this way, 
the worst possible error in computation for the metric 
potentials of BH is estimated. Even near the radius for 
a = 0, the second order convergence is maintained as 
seen in the figure. We also show in the bottom panel of 
Fig. [3l averaged fractional errors for the conformal fac- 
tor if). The value of tp is about 2 near r = r a , and for 
such a potential the convergence of the solution is almost 
uniform in all radii as observed. 
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FIG. 4. Same as Fig. [3] but for the equal mass BBH data 
calculated on a single patch Fig. [5] Top panel: Plots for the 
fractional errors in the lapse 8a /a for the equal mass BBH 
data along the positive s-axis. Middle panel: averaged frac- 
tional errors in the lapse (\8a/a\). Bottom panel: averaged 
fractional errors in the conformal factor (\8ip/^\). Plots show 
the errors with changing the number of all grid points as Dl- 
D4 in Table [ml 





FIG. 5. A setup for a single coordinate grid patches for a cal- 
culation of equal mass BBH. The radius of coordinate patch 
doesn't reflect the actual size. 



Type Patch 


r a 


n 


r c 


r e 


iV r f iV, m 


iV r 


N e 




L 


El 


COCP-1 


0.2 


10 4 


1.25 


1.125 


16 


20 


64 


24 


24 


12 




COCP-2 0.4 


10 4 


1.25 


1.125 


16 


20 


64 


24 


24 


12 


E2 


COCP-1 


0.2 


10 4 


1.25 


1.125 


32 


40 


128 


48 


48 


12 




COCP-2 0.4 


10 4 


1.25 


1.125 


32 


40 


128 


48 


48 


12 


E3 


COCP-1 


0.2 


10 4 


1.25 


1.125 


64 


80 


256 


96 


96 


12 




COCP-2 


0.4 


10 4 


1.25 


1.125 


64 


80 


256 


96 


96 


12 


E4 


COCP-1 


0.2 


10 4 


1.25 


1.125 


128 


160 


512 


192 


192 


12 




COCP-2 0.4 


10 4 


1.25 


1.125 


128 


160 


512 


192 


192 


12 


Fl 


COCP-1 


0.2 


10 2 


1.25 


1.125 


16 


20 


48 


24 


24 


12 




COCP-2 


0.4 


10 2 


1.25 


1.125 


16 


20 


48 


24 


24 


12 




ARCP 


5.0 


10 6 


6.25 




4 


5 


48 


24 


24 


12 


F2 


COCP-1 


0.2 


10 2 


1.25 


1.125 


32 


40 


96 


48 


48 


12 




COCP-2 


0.4 


10 2 


1.25 


1.125 


32 


40 


96 


48 


48 


12 




ARCP 


5.0 


10 6 


6.25 




8 


10 


96 


48 


48 


12 


F3 


COCP-1 


0.2 


10 2 


1.25 


1.125 


64 


80 


192 


96 


96 


12 




COCP-2 


0.4 


10 2 


1.25 


1.125 


64 


80 


192 


96 


96 


12 




ARCP 


5.0 


10 6 


6.25 




16 


20 


192 


96 


96 


12 


F4 


COCP-1 


0.2 


10 2 


1.25 


1.125 


128 


160 


384 


192 


192 


12 




COCP-2 


0.4 


10 2 


1.25 


1.125 


128 


160 


384 


192 


192 


12 




ARCP 


5.0 


10 6 


6.25 




32 


40 


384 


192 


192 


12 



TABLE IV. Grid parameters used in convergence tests for 
non-equal mass BBH data solved on multiple coordinate 
patches. The separation of two BHs is set as d s = 2.5. 



3. Non equal mass BBH computed with multiple patches 



Convergence test for non-equal mass BBH data using 
multi coordinate patches discussed in Sec|H]is performed 
with grid parameters presented in Table IIV1 In the first 
example. BBH data is computed on two COCP whose 
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r 

FIG. 6. Same as Fig. [3] but averaged fractional errors in the 
lapse (|<5a/a|) are plotted along the radial coordinate r for 
the non-equal mass BBH data calculated on multiple patches 
Fig. [1] Top panel: data computed on two multiple patches 
with changing the number of grid points as E1-E4 in Table 
Mil Bottom panel: data computed on three multiple patches 
with changing the number of all grid points as F1-F4 in Table 

EU 



boundary radius r = r& is taken large enough to reach to 
asymptotics = 10 4 . In this computation, the number 
of grid points is chosen as E1-E4 in Table IIV1 in which 
the resolution in radial grids in the region r > r c is in- 
creased by 44/24 times of the corresponding level of res- 
olution for equal mass BBH case, D1-D4. Separation be- 
tween the coordinate centers of two BH is set as d s = 2.5, 
the excision radius and mass parameter of the first BH 
are r a = 0.2 with Mi = 2r a = 0.4, and those of second 
BH r a = 0.4 with M 2 = 2r a = 0.8. The results of the 
averaged fractional error in the top panel of Fig. [5] is sim- 
ilar to those of the equal mass BBH case in Fig. @J which 
proves our multiple patch methods works accurately as 
expected. 

Finally, the BBH data is computed using three multi- 
ple patches shown in Fig. [TJ The number of grid points 
is chosen as F1-F4 in Table HVl and the separation is set 
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FIG. 7. Plots for the conformal factor tp and the lapse a, 
computed on the three multi patches. The model is the same 
as that of the bottom panel of Fig. [5] 
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o s 


Spin axis 


Figures 


TU 


3.0 


1.0 


0.3 


0.0 




Fig. [9] 


TU 


2.2 


0.005 


0.08 


0.0 




Fig. [10] 


AH 




0.1 


0.08 


0.0 




Figs. [T2lfT8l 
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0.08 
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0.1 


0.08 


0.1 


z-axis 


Fig. El 


AH 




0.1 


0.08 


0.1 


y-axis 


Fig.UJ 



TABLE V. Boundary conditions and their parameters used 
in the computations for BBH initial data. The first column, 
Type, denotes the types of boundary conditions used, TU 
corresponds to Eq. (|35p . and AH corresponds to apparent 
horizon boundary conditions Eqs. (|36I) - (I38I ). The model of 
Fig|9]is computed with a binary separation d s = 2.8, a radius 
of BH excision r a =0.1, and a radius of binary excision r e = 
1.3. All the other models are computed with the parameter 
set D3 in Table ITTTI 



as d s — 2.5. In the bottom panel of Fig. [5J the results for 
the averaged fractional errors (\6a/a\) are shown. In this 
computation we decreased the values of mass parameter 
as Mi = 0.8 x 2r a = 0.32 with r a = 0.2 for the first BH 
and M 2 = 0.8 x 2r a = 0.64 with r a = 0.4 for the second 
BH, so that the lapse a is always positive even near the 
BH excision boundary at r — r a . As seen in the figure, 
the error near the BH boundary is decreased about 1/10 
of the previous case, although the resolutions near the 
BH are the same. In the Fig. [71 we present the plots 
of the potentials a and ip for the same model to show 
a smooth transition of potentials from one to the other 
patch. 
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FIG. 8. Schematic figure for the n rotation symmetry of x and y components of the shift from a region inside to S' e to the 
excised region inside of the sphere S e - 




FIG. 9. Plots for the y-component of shift vector along x-axis 
of BBH initial data computed from the boundary conditions 
(|35[) . Parameters in the conditions are chosen the same as a 
solution presented in [3 as m = 3, no = 0.1, and SI = 0.3. 
Two BH are located at x = and x = 2.8. The region 
inside of the excised sphere S e is interpolated using 7r-rotation 
symmetry. Thin solid (black) vertical lines are the boundaries 
at S a with the radius r a =0.1, and dotted (blue) vertical lines 
are the boundaries at S e with the radius r e = 1.3. 



FIG. 10. Plots for the conformal factor ip and lapse func- 
tion a along the a;-axis of BBH initial data computed for the 
boundary condition (|35[) . Parameters in the conditions are 
chosen as ni = 0.1, n = 0.005 and SI = 0.08. The two BH 
are located at x — and x — 2.5. The region inside of the 
excised sphere S e is interpolated using 7r-rotation symmetry. 
Thin solid (black) vertical lines are the boundaries at S a with 
the radius r a = 0.2, and dotted (blue) vertical lines are the 
boundaries at S e with the radius r e = 1.125. 



IV. INITIAL DATA FOR BINARY BLACK 
HOLES 

Finally we present examples for initial data sets for 
equal mass BBH, which have been widely used for BBH 
merger simulations in the literature. Among several for- 



mulations for computing such data sets (see e.g. [H@ and 
references therein), we adopt Isenberg-Wilson-Mathews 
(IWM) formulation. In this section, we show the solu- 
tions computed from two different types of boundary con- 
ditions. The first are simple boundary conditions used 
in our previous paper [l4j . The second are the appar- 



14 



eo. 



0.03 
0.02 



-0.01 
-0.02 
-0.03 



0.15 
0.10 
0.05 
odT 0.00 
-0.05 
-0.10 
-0.15 

0.15 
0.10 



- 


' COCP BH 




i 








> 




— 














- 




















r 


















r- ■ 






























































l 


, 







0.05 



0.00 





y 







i 










' 1 , 


COCP-BH 




































































































l 















-6 


-4 


-2 





2 


4 


6 






■ 


< 




i 

COCP-BH 




i 




i 


i 


































































l 




, 


, 







y 



FIG. 11. Plots for the same model as Fig. [10] but for the 
components of the shift j3i. Top panel: fi x along the y-axis. 
Middle panel: fi y across the rr-axis. Bottom panel: /3 y across 
the y-axis. 



ent horizon boundary conditions, which have been used 
to compute quasi-circular initial data for BBH (see e.g. 

B B HE 13). 
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FIG. 12. Shift vector on the xy-plane of the BBH initial data 
for the case with AH boundary conditions (|36[) - (|38[ ). Param- 
eters in the conditions are chosen as no = 0.1, Q = 0.08, and 
fi s = 0. The center of mass is located at x — 1.25 on the 
a;-axis. In the computation, the region inside of the left thick 
circle centered at the origin with a radius r a = 0.2, and of 
the thin circle (green) with a radius r e — 1.125 centered at 
x = 2.5 on the a>axis are excised. The data inside of the thin 
circle is interpolated by a symmetry. Note that the center of 
mass does not coincide with the center of the x-axis. 



A. IWM formulation and boundary conditions 

IWM formulation has been widely used for construct- 
ing quasi-cquilibrium initial data for binary compact ob- 
jects. We summarize the basic equations below (for more 
details see e.g. @, [H, H3])- The spacetime metric on S t 
is written in 3+1 form as 

ds 2 = g lll/ dx p, dx v 

= -a 2 dt 2 + 7^ {dx i + [3 l dt){dx ] + /3 j dt), (30) 

where the spatial three metric 7^- on the slice S t is as- 
sumed to be conformably flat 7^ = ip fij. Here, field 
variables tp, a, and j3 l are the conformal factor, lapse, 
and shift vector, respectively. We also assume maximal 
slicing to St , so that the trace of the extrinsic curvature 
Kij := — 5^(£*7ij — ^filij) vanishes. Then, writing its 
tracefree part Aij , the conformally rcscalcd quantity Ay- 
becomes 



1 

2a" 



(31) 



where the derivative di is associated with the flat metric 
fij, and conformally rescaled quantities with tilde are 
defined by = A^ and fi l — whose indexes are 
lowered (raised) by fij (f 1 ^)- The system to be solved, 
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FIG. 13. Contour plots on the xy-plane for the conformal 
factor tjj (top panel) and the lapse function a (bottom panel) 
for the same model as Fig. 1121 For if), we draw isolines from 
ip = 1.1 to "0 = 2.0 with step 0.1. For a, we draw isolines 
from q = 0.2 to ot = 0.9 with step 0.05. 



which arc Hamiltonian and momentum constraints and 
the spatial trace of the Einstein's equation, becomes 

v> 5 ~. 



-AijA 3 , 



(32) 



A/?, = -2 a A^d, In ^ - V , (33) 

a 3 



A(cn/>) = -a^AijAv, 



(34) 



where A := did 1 is a flat Laplacian. It is noted that, for 
the shift equation Q33p. the Cartesian components are 



solved on the spherical coordinates. 

As a first set of the boundary conditions at the BH 
excision boundary S a and at the boundary of computa- 
tional domain Sb for the above system Eqs. (|3"2"|) - (f3"4"|) . we 
choose the following for simplicity, 



V' I T—T 



"1 



and 



n 



r=r b 
\r=r b 
\r=r b 



= 1.0 

= 0.0 (35) 
= 1.0 



where n\ and no are arbitrary positive constants taken as 
ni > 2 and no < 1, 4> l cm = (-J/™, x cm , 0) is the rotational 
vector with respect to the center of mass of the binary as- 
sociated with coordinates (x cm , y cm , z cm ) := (x — d, y, z), 
and corresponds to the orbital angular velocity. The 
radius r\, is taken large enough as in the test problems 
in Sec. IIIII Despite the fact that the boundary condi- 
tions above may be of the simplest type for IWM for- 
mulation deduced for acquiring BBH data with non-zero 
orbital angular momentum in an asymptotically flat sys- 
tem, they capture the qualitative functional behavior of 
the unknown fields {ip, a, /?*}, as more realistic boundary 
conditions mentioned below. The solutions calculated 
from these boundary conditions (|35|) are compared with 
our previous code [14[ that uses a different structure for 
the multiple spherical coordinate patches, as well as those 
solutions of different boundary conditions. 

For a second set of the boundary conditions, we impose 
more realistic boundary conditions at the BH boundary 
S a , in particular those that represent apparent horizons 
in equilibrium @, i, H, S3 , 



dip ip 
dr 2r 



4 

n i n li 
= "0, 



(36) 

(37) 
(38) 



where no is an arbitrary positive constant for which we 
choose no < 0.1, s l is the unit normal to the sphere S a , 
and il s represents the spin of each black hole. The vector 
<j>l is the rotational vectors with respect to the coordi- 
nate center of BH that generates the BH spin. The spin 
axis is not necessary parallel to the z-axis. Demanding 
the sphere S a to be an apparent horizon (AH) results in 
Eq. p6p . while demanding the horizon to be in equilib- 
rium results in Eq. (1 3 7 1) . 

For the present calculations for BBH initial data, we 
also assume 7r-rotation symmetry of the system around 
the center of mass. That is, two BH have equal masses, 
and 7r-rotation symmetric spins if any. In other words, 
the same boundary conditions are imposed on both 
BH. In those cases, a single patch method discussed in 
Sec. lIIIB~T1 can be used for simplicity. As shown in PiglHl 
the metric potentials are mapped to the excised sphere 
S e from the corresponding sphere S' e , taking into account 
the parity of the variables with respect to the 7r-rotation. 
As an example, the 7r-rotation symmetries of the shift 
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components on the xy-plane are shown schematically (the 
shift at A' is mapped to point A), together with the cor- 
responding rules for the derivatives of a function along 
the x-axis inside the excised sphere (B'C mapped to 
BC) and along the y-axis (B'D' mapped to BD). In 
terms of the center of mass coordinates, the 7r-rotation 
symmetries of the components of the shift vector become 

fix ( *^cni5 Vcim ^cm) ^(^cmi Vcrai ^cm) ) 

$y\ ^cmi Ucm: ^cm) ftyi^cTm ?/cm> ^cm)- 

The z-component j3 z is mapped as a scalar quantity. 

Mapped quantities are used in the elliptic equation 
solver when the sources of surface integrals on the ex- 
cised sphere S e , Eq. ([7]) are evaluated. Also, at the end 
of each iteration step, obtained potentials {ip, a, fa} be- 
tween the spheres S a and S' e are interpolated inside the 
sphere S e following the same rules for the parity of each 
variable. At the end we have the solution at every point 
inside St, and outside of the two black holes of radius r a 
positioned at the origin and at x = d s on the x-axis. 

In the case of no spins f2 s = 0, or spins are parallel 
to one of coordinate axis, additional symmetries with re- 
spect to the xy\ cra plane occur. In our previous codes 
[12 - [l4| . a part or all of these symmetries were encoded 
in the elliptic solver, for a computational domain was re- 
duced by assuming the symmetries so that only a part of 
the whole hypersurface T, t was solved. In COCAL, we are 
solving in the whole E ( and therefore such symmetries 
are satisfied in the solution within negligible numerical 
errors. 



B. Solutions for BBH initial data 

Finally, we present the BBH initial data sets computed 
from the above formulation with several parameter sets 
for boundary conditions listed in Table |V] As we con- 
centrate on testing the COCAL code, we do not discuss 
much of the physical contents of the initial data sets, but 
display the plots of the fields to check their behaviors. 

In Fig. [51 we calculated BBH data for the same model 
as shown in our previous paper [l4| for a comparison. 
Parameters in the boundary condition ([35)1 are chosen as 
n\ = 3.0, no = 1.0, and £1 = 0.3. Only in this model, 
we choose the binary separation as d s = 2.8 and the BH 
excision radius r a = 0.1. We find the solution agrees 
well with Fig. 11 of [l4[ as expected, although the struc- 
tures of coordinate patches are different in each code. In 
Figs.[Tni and [Til We use the same boundary conditions 
(|35[) but with different parameters which may be more 
common values for the BBH data. In the computation, 
the grid parameters used are D3 in Table IIIII When the 
value of f2 is increased, the magnitude of the lapse and 
the conformal factor remain almost the same, while the 
magnitude of the components of the shift increase with 
the functional behavior staying the same. For example 
when fl = 0.1, f3 x on the y-axis varies between ±0.05 
while j3 y on the x-axis varies between ±0.15. For this 



boundary conditions, the code blows up when fl is ap- 
proximately greater than 0.2. 

Solutions for the BBH initial data with the AH bound- 
ary conditions (|3"6")) - (j3"5)) are shown in Figs. H^HTBI In 
this calculation, the resolution is D3 in Table IIIII with 
il = 0.08 and il s = 0. The shift vector plots and the 
contour plots of the conformal factor ip and the lapse a 
in xy-plane are for the model with uq = 0.1. The behav- 
ior of the field ip and a are analogous to those from the 
simple boundary condition ([33]) shown in Fig. [TO] Be- 
cause of the choice fi s = 0, the solution satisfies xy-plane 
symmetry. 

From a comparison between the results of parameters 
?io = 0.1 and 0.005 shown in Figs. fTHfTBl it is found that 
the results with no = 0.005 becomes more similar to the 
results of the first boundary conditions (|35[) shown in 
Fig llll For example, it is most evident in the plot for /3 y 
along the y-axis, FigfTS] middle panel and FigfTTI bottom 
panel. This seems to be a correct behavior because the 
first term of Eq. (|37j) contributes less as the value of alpha 
(parameter no) become smaller 0. However, it is not 
expected that there is a solution in the limit of no — > 0, 
in fact, for both types of boundary conditions iterations 
diverge when the no < 0.004 in the COCAL code. 

All of the solutions above have zero black hole spins. 
Setting spins f2 s = 0.1 in the same direction as the or- 
bital motion (rotation on the xy-plane) we get a solution 
similar to Figs. [T4TfTBl except for f3 x along the y-axis and 
/3 y along the x-axis. The results with and without spins 
are compared in Fig. [T7] for the case with no = 0.1. For 
more general black hole spins we obtain Fig[TS] where we 
have taken a spin Q s = 0.1 along the y-axis. In these 
plots for the solutions with the spins, we confirm that all 
components of the shift vectors behave correctly along 
each coordinate axis. 



V. DISCUSSION 

Although the numerical method presented in this pa- 
per may seem similar to the one presented in our previous 
paper [lj| , the robustness of the convergence and the con- 
trol of the numerical errors are largely improved. In some 
cases, the previous method failed to compute a contin- 
uous solution at the interface between multiple patches 
during the iterations, and therefore a convergence to a 
solution couldn't be achieved. A major reason for this 
failure turned out to be a lack of enough overlap re- 
gion between coordinate patches. The numerical errors of 
the field variables are relatively larger near the boundary 
of computational domain as seen, for example, in Fig|5] 
Hence if the overlap region is small, those potentials with 
larger numerical errors overlap, which seems to cause the 



4 For the first boundary conditions II35H . the solution of the shift 
does not depend much on no in the region no < O.i. 
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non-convergence of the iteration. In the COGAL code, 
the overlap region is almost as large as the whole domain 
for the two coordinate patch configuration, and is large 
enough even for the three coordinate patch configuration. 
We have never observed so far a discontinuous behavior 
of a field in the solutions of COCAL. 

The COCAL code currently runs only on a serial proces- 
sor, which is sufficient to maintain the accuracy presented 
in Sec lIIII In the computation for BBH initial data shown 
in Sec. IV, the size of main memory and CPU times per 
1 iteration cycle used by COCAL are about 800MB 50sec 
for D3 grid, and 6GB 8min for D4 grid, and around 50- 
150 iterations are needed for a convergence, where the 
iterations start from an initial guess ip = a = 1 and 
Pi =0. Because we use second order accurate formulas 
in COCAL, we can decrease the numerical error by two 
orders of magnitude with 10 times more grid points in 
each direction, that is, 10 more grids in total. Con- 
sidering specs of common parallel computer systems it 
seems to be feasible to achieve this accuracy by paral- 
lelizing COCAL. We have started to develop a prototype 
of such parallelized COCAL code, whose results would be 
presented elsewhere. 

The most advantageous feature of COCAL would be its 
simplicity in coding. This helps the users to introduce 
more complex physics on top of the current code. For 
example, we have developed subroutines for solving spa- 
tially conformally flat data (IWM formulation) first, and 
later added subroutines for solving non-conformally flat 
data (waveless formulation) on top. In the same way, 
it will be straight forward to incorporate subroutines to 
solve electromagnetic fields, which enables us to investi- 
gate for example magnetar models. We proceed to de- 
velop codes for computing various kinds of astrophys- 
ically realistic equilibriums and quasi-equilibrium data 
and provide the results to applications including initial 
data for numerical relativity simulations. We plan for 
making the COCAL code and computed initial data sets 
available for public in the near future. 
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Appendix A: Computations for homogeneous 
solutions 



In this Appendix, we show a concrete derivation for 
the homogeneous solution xi x ) used in the elliptic solver 
([p]) for two cases, one with Neumann boundary condition 
at the inner boundary sphere S a and Dirichlct bound- 



ary condition at the outer boundary sphere Sb, and the 
other with Dirichlet conditions at both S a and Sb- Wc 
summarize the other cases in the next Appendix [Bl 

As explained in Sec. Ill C[ the solution of = S is 
written $(:r) = x( x ) + ^int(^) to impose certain bound- 
ary conditions at the two spheres S a and Sb- In order 
to do so, the homogeneous solution x of the Laplacian 
is split into two functions Xa{x) and Xb(x) as x( x ) — 
Xa(x) + Xb(x). Both are solutions of Laplace equations, 
and one for the exterior of the sphere S a and the other for 
the interior of the sphere Sh (see Fig. [19)1 . Since r l , r~ l ~ x 
are the solutions of the radial part of the Laplacian, the 



contribution Xa is taken to be a series of \ 



while the 



contribution Xb to be a series of r . Therefore, wc write 



oo £ 



Xa{x) 



1 \ - \ - (£-m)\ . . _ 
4^ e ^ (c ° S0)r 

1=0 rn=0 v ' 

x [Ae m cos(m</>) + B im sin(m^)] 



1 ^ (£-m)\ 
'(i + m)\ 



Pf (cos 9)r £ x 



47T 

1=0 m=0 

x \p£ m cos(m</>) + Di m sin(m</>)] 



(Al) 



(A2) 



where Ai m , Bi m , Cg m , and Dg m arc constants. One com- 
mon choice is 



<9$ 

dr 



BC 



dr 



(A3) 
(A4) 



where \ r =r a and $Bc|r=ri, are known functions de- 
fined at boundaries S a and Sb, respectively; we consider 
the case with Neumann boundary condition at the inner 
surface S a and Dirichlet at the outer one Sb- 

From boundary condition Eq. (|A3)I with the use of 
Eqs. (jATj) . (fA2|) . @, 0, and the orthogonality relations 

2 (f + mV 
i Pr(co S 6)PF(co S 6) S med6 = ^ + 1 

2-7T 

sin(m<^) cos(m' 4>)d(j) = 

o 



cos(to0) cos(m' cf>)d(f> 



2tt 



we get 



,.''+2 



■A im +£ri- 1 C tm = (2e+l) x 



7T />2?r 

o Jo 



9$ 



BO 



INT 



dr 



dr 



Pt{cos9) cos(m(t))dn 



(A5) 



1 



-Be 



r, 

7T p27T 



1+2 ^ a 

1 

( <9$BC 9$INT 



t _! 2(21+1) 
r„ D em = — x 



oh V dr 



dr 



PP(cos0)sin(m(j>)dn 
(A6) 



18 



Using now boundary condition Eq. 



4J) and again the same equations we get 

r^- l A im + riC lm = {21+1) x 

n27T 
($bc - $iN T ) r=rb PP(cos6) cosMdn (A7) 

„ 2(2* +1) 



n27T 
($BC - * 



iN T ) r=r Pr(cos#)sin(m</>)dn (A8) 



The system of four equations with respect to Ai m , Bg m , 
Cim, Dim can be solved to yield 



.4 



tin 



a 



D 



i III 



i+1 \r b 



2£+l 



e(2£+l) v+1 fr a 



£+1 a 

e + i a 



—\ I ($bc - $iNT)^ m (cos6»)cos(m^)^ 

r bJ J Sb 



<9$bc d$ 



INT 



dr dr 



PP(cos9) cos(m0)drj, 



£+1 \r b 



2£+l 



{2£+l)b- e [ ($bc - *iNT) J P" i (cos6»)cos(m ( /))drj 



2i+l r 2 a fr a \ e f /a$ BC 9$ 



t + 1 rf. +1 \ r b J J s V dr 



J INT 



dr 



PP(cos9) cos(m(/))dn, 



' + 1 \n 



2£+l 



2£(2e+l) t+1 (r a 



{£+l)e m a \r b 



($bc - $iNT) J P f m (cos6») sm(m(j))dn 



s b 



t + 1 \n 



21+1 



(e + i)e m a 

2{2£ + l) _ t 



+ 



e m Js b 

2{2t + l) r\ (r a 



dr dr 
($bc - $int)-P™(cos0) sin(m0)dfi 



0? + l)e m r£ +1 W 
Substituting the above in Eqs. (|A1[) and (|A2[) . we have 
1 j (9 — \\ 

1=0 m=0 V 7 ' 



/<9$BC <9$ 



INT 



dr 



PP(cos6) sin(m0)cftl 



(21+1) - 

n 



t+1 \ r J 



1 + la 



2£+l 



+ - 



-2^-1 



>M-1 



'+1 Vn, 



2£+l 



[*bc - $iNT]iT(c os6,/ ) cos[m(0 - 

<9$BC 9$INT " 



9r 9r 



P f m (cos6»')cos[m(0-(/)')]^' > . (A9) 



r 



Similarly, when we have Dirichlct boundary conditions the contribution to the potential will be 
on both the inner and the outer spheres S a and 5b, 



^*lr— r a ^*Bc| r — Ta 1 
^lr=n, — ^BC | r=ri) j 



(A10) 
(AH) 
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Xa(x) +Xb(x) = 



^ oo £ 

E E 



47T ^— ' (« + ml! 

£=0 m=0 v y 



(2£+l) 



, 2£+l 



-{21+1) — 



£+1 (Ik. 



-1 



2£+l 



[$bc - $int]JT(cos6I') cos[m(</> - </>')]dn' 



[$bc - $iNT]P f m (cos 9') cos[m{4> - $')]dSl' 



> (A12) 



The final solution will be obtained from the iteration of 



= X (a?) + $int(.t) = Xa(x) + X b(x) + $int(z) Xb = ^ (E [fff'forOM - 3 r / 5 ? C (r, r') S 

"'it £=0 



and 

Xb 



x 

r'=r b 



where \ a + \b are taken from Eq. (|A9|) or (|A12[) depend- 
ing on the boundary condition. For the other boundary 
value problem, x( x ) = Xa{x) + Xb{x) is modified accord- 
ingly as shown in the next Appendix. 



Appendix B: Green's functions and surface integrals 

In this Appendix, we present the explicit forms of the 
kernel functions denoted by G {x, x') in Eq. (fl~2j) . which 
appear in the surface integrals of the homogeneous so- 
lution x( x ) CD- Various types of boundary conditions 
are imposed on a spherical domain bounded by two con- 
centric spheres S a and Sb, and the corresponding kernel 
functions available in COCAL are tabulated in Table Q] In 
[3 we used Green's functions G NB {x,x'), G DD {x,x'), 
and G ND (a;, x 1 ). Here we construct also G DN {x, x'), 
G NN (x, a;'), G RD (x,x'). 

The surface integral 



x{x) = h 



S a US b 



G BC {x,x , )V' a ^{x') 
${x')V' a G BC {x,x')] dS' a , (Bl) 



where <&(V) := $ BC (x') — $int(^'); are written below for 
both the exterior and the interior problem. Noticing that 
dS' a is pointing outward, S a and Sb are the concentric 
spheres, and V a f{x') dS' a = d r >fr /2 dQ', we have 



Xa — 



4-7T 



r. OO 

/ Yu [-9t C (r,r')d r ^ + d r ,gf c (r,r')$ 



E 



- m V. ,1 Pt n (cos6)Pl n (cos6')coslm(4>-(!)')]r 2 a dn', 



(e + m)\ 



(B2) 



E 



(£-m)\ 
'{l + m)\ 



PP(cos6)Pp{cos0')cos[m((f> - (j)')]r 2 b dQ,' . 

(B3) 



Note that these Xa and Xb in this section are defined dif- 
ferently from those in the previous section, Eq. (|A1|) and 
(|A2j) . As shown in Fig. \T9\ we denote by r a the radius 
of the sphere S a for the exterior problem and by Xa the 
corresponding integral and by r& the radius of the sphere 
Sb for the interior problem and Xb the corresponding in- 
tegral. 



1. Kernel function G NB {x,x') 

The radial part for the kernel function without bound- 
ary G NB (x,x'), is 



NB/ /\ _ ' < 



\r,r') 



(B4) 



where r> := max{r, r'}, and r < := min{r, r'}. The radial 
part of the kernel function in the surface integral on S a 
(IB2D becomes 



9i(r,r a ) = «>,r a ) = (B5) 



while the surface integral on Sb 



9£ B (r,r b ) = — , d r ,g? B (r,r a ) = -{1+1)- 



r b r b 



(B6) 



2. Kernel function G DD {x,x) 



When the Dirichlet boundary condition is imposed on 
both S a and Sb, the kernel function G DD (a;,a;') satisfies 



G uu {x,x')\ Sa =0, G DD (x,x')|, 
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These conditions lead to the vanishing of the radial part 
gf B {r, r') on the two spheres S a and Sb, 

9? D (r,r a )=g? D (r,r b ) = 0, 

which result to following formula for gf D (r, r'): 



DD/ f\ 

9i [r, r ) 



1 - I — 

n 



21+1 



e+i 



l+l 



n \ ( r> 



The radial kernel function at S a becomes 
gf D (r,r a )=0, 



.(B7) 



/ e+i / x2£+i 



Of) 



r-. (B8) 



and at Sb, 

gf B (r,r b ) = 0, 

d r .gf D (r,r b ) = -(2£+l) 



3 ( -( B9 ) 



J!+2 



Special cases are when the surface Sb is absent in the 
limit Tb — > oo, 



d r ,g? D (r,r a ) = (21+1)^- 



(BIO) 



for the surface integral at S a , or when the surface S a is 
absent in the limit r a — > 0, 



ft-'flf (r,r t ) = -(2* + l)- 



.«+2 : 



(Bll) 



for the surface integral at Sb- The latter will be used for 
computing neutron stars. 



3. Kernel function G NU (x,x') 

When the Neumann and Dirichlct boundary conditions 
are imposed on S a and Sb, respectively, the kernel func- 
tion G ND (x,x') satisfies 



d r >G ND (x,x')\ Sa = 0, G ND {x,x')\ Sb 
or in terms of gf B (r, r') 







Then the radial part of G (x,x') is 



£ 




l+l 


- 


\ 1 




t) 


+ £ 


-t 


-l 


r >j 





2£+l 



J+ 1 



n 

The radial kernel function at S a becomes 



(B12) 



2£ + l ri VT 



£ + 1 r' +1 



21+ 



r ,(B13) 



6 1 + _L_ (Z. 

1 ^ e+i \r bj/ 

dr'gf D (r,r a )=0 
and at S b , 

gf D (r,r b ) = 



^flf u (r,r6) = -(2£+l) 



..'+2 



1 



f+1 (r-b) 



2£+l 



(B14) 



Special case is when the surface S b is absent in the 
limit rb — > oo, 



2^+1 



(B15) 



for the surface integral at S a . 

4. Kernel function G DN (a:,2:') 

When the Dirichlet and Neumann boundary conditions 
are imposed on S a and Sb, respectively, the kernel func- 
tion G DN (a;,a; / ) satisfies 



G™(x,x')\ Sa =0, d r ,G DN (x,x')\ Sb 
or in terms of g^ N (r, r') 

g? N (r,r a ) = d r ,g? N (r,r b ) = . 
Then the radial part of G BN (x,x r ) is 

t 







2f+l 



f+1 



£+1 \r> 



e+i 



(B16) 
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The radial kernel function at S a becomes 
<?? N (r,r o ) = 0, 

9r , 9 ?»(r,r a ) = +t+A :.l 



1+ 



2e+i > 



(B17) 



and at Sb, 



£+1 r{ +1 _L 
£+ 



^fl? N (r,r t ) = 0. 



5. Kernel function G RD (a:,2;') 

When the Robin and Dirichlet boundary conditions are 
imposed on S a and Sb, respectively, the kernel function 
G RD {x,x') satisfies 



g G RD G RD 



dr 2r j g 

or in terms of <7^(r, r') 



0. 



G RD (x,x') 







0. 



5 f D (r,r b )=0. 



9r 2r 
Then radial part of G RB (x, x') is 



9? U (r,r') 



e+i 



21+1 



J+ 1 



.(B19) 



For the surface integral at S a , it is more convenient to 
rewrite Eq.(|B2]). 



Xa = h L £ 



_,RD 



E 

m= 



+ (^# D (r,r') + ^^ 



(7 - mV 

Wt: J -PP(cos6)Pp(cos9') cos[m(0 - ^Ol^dfi', 



(£ + m)! 



(B20) 



while for Xb, Eq. (|B3|) is used. Here, the radial kernel 
function at 5 a becomes 



M-i 



1+^ 



2e+i 



(B21) 



^fl< (r,r a ) + 



gf u (r,r a ) 

2r a 



The radial kernel function at Sb in Eq. (|B3[) becomes 

gf D (r,r b ) = 0, 



d r ,g? u (r,rb) = -(2£ + l)-£ 



a+i 



1+2 
'b 1 



(») 



Special case is when the surface Sb is absent in the 
limit rb — > oo. In that case, we use in Eq. (|B20j) . for 
£=1,2,... 



flT(r,r.) = 2^. 



6. Kernel function G DK {x, x') 



(B23) 



When the Dirichlet and Robin boundary conditions are 
imposed on S a and Sb, respectively, the kernel function 
G DR (x,x') satisfies 



G DR (x,x')\ 



0, 



9G DR G DR 



9r 



2r 







or in terms of gf R (r, r') 



g? K (r,r a ) = 0, 



d9? R , 9? R 
dr 2r 



Then radial part of G BK (x,x') is 



9? R (r,r>) 



[ - 

n 



2i+i 



e+i 



.(B24) 



For the surface integral at Sb, it is more convenient to 
rewrite Eq. (|B3|) , 



Xb 



DR/ I 



<I> 



(r,r')(^$ + — 



(£-m)\ 



2r' 



x 

r'=r b 



V e m ^ "^> f (cos g)i^ (cob ^) cos[m(^ - </>')}r 2 b dn' , 

(B25) 

while for %„, Eq. (|B2|) is used. Here, the radial kernel 
function at S a becomes 



g^{r,r a ) = 0, 



0. 



9 r ^ K (^^) = (2^ + l)i+T 

r <> 1 



(r.y +(tlL) W 
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The radial kernel function at S b in Eq. (|B3[) becomes 

2— (B27) 



t (A 



d r ,g? R (r,r b ) 



'b 1 



9r(r,r b ) 
2r b 



0. 



Special case is when the surface S a is absent in the 
limit r a — > 0. In that case, we use in Eq. (|B27|) . for 
1 = 1,2,... 



_,DR 



(r,r b ) = 



(B28) 



Such kernel functions for imposing Robin boundary 
conditions at the outer surface S b may improve an ac- 
curacy of the solution especially near the boundary [30j . 



7. Kernel function G NN {x,x') 

When the Neumann boundary condition is imposed on 
both S a and Sb, the kernel function G NN (a;, x 1 ) satisfies 



c^G NN (x,x')| 



G a , 



d r ,G NN (x,x')\ 



Gh 



where G a , Gb cannot be both zero. In that case 
G NN (x,x') does not exist since the £ = mode cannot 
be satisfied. Therefore the boundary conditions for the 
radial part will be 



dr'9e {r,r a ) = G a 5 oe , d r >g e (r,r b ) = G b S M 
Then for £ = 1,2,... we get 

2£+l~ 



1 



£ + 1 



£ + 1 

i 



J+ 1 



£+i y r> 

(symmetric in r,r'), while for £ = 

G r 2 



?+l 



(B29) 



1 



r> r 

where h(r') arbitrary function. Symmetry is imposed by 
choosing h(r') = —G a r 2 /r' therefore 



9 r(ry) = --riG a (± + ± 



(B30) 



Note also that in order to satisfy the £ = mode the 
following condition must hold 

1 = G a r 2 — GbT 2 



therefore G a and Gb cannot be chosen arbitrarily. The 
surface integral at S a is 



oo t 



^ = iEE4^^r(cos, 



4tt ^ '""(l + mY.' 



J (-9f N (r,r a )^ + G a 6 oe ® \ x 
xP7 l (cos6»') cos[m{((> - <f>')}iidn' (B31) 



where for £ = 1,2, . 
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_1_ t (r b \£+l 


2£+l r e a 


(*) 


i+1 \ r J 



£ rf +1 



2e+i 



(B32) 



and £ = 



3r(r 1 r ) = --rX( 1 + 1 



The surface integral at Sb is 

1 (£-m)\ 



(B33) 



4tt ^ ^ + m)!" 



P7 1 (cos 6) x 



x^ (g^ {r ,r b )^--G b S ^ x 

xPp(cos0')co8[m(<f> - 4>')]rldQ! (B34) 
where for £ = 1,2,... 
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£ rf +1 



1- ^ 



and £ = 



.9o NN (r,^) = 1 -^G a (I + i 



(B36) 



Special case is when the surface 5 a is absent in the 
limit r a — > 0. In that case, we use in Eq. (|B34[) . for 



1,2, 



and £ = 



2^+1 



9o NN (^r 6 ) 



(B37) 



(B38) 



Appendix C: Finite difference formulas 

The second order finite difference formulas used in the 
elliptic equation solvers of COCAL code are summarized 
in this section. In evaluating the integrals of the solver, 
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such as Eq. ([7]), we use the mid-point rule. Hence, we 

need to evaluate the source terms in the integrand at the d(j> l_ 5' J _ 2' 2 



mid-points (r 4 _i,^_i,0 fc _i) of the grid points that may 1 i j ffaBj^) - fin^fa-t) 

involve values of potentials and their derivatives. Those — — 

are calculated, respectively, by r=i-i J=j-i 

- 5 E E E frTMK), (CI) 



E E 

Z=i-1 J=j-1 if=fc-l 

-J-(ri_i.,9 4 _i,6 k _i) 
Or 2 3 2 2 



The quadrature formula for the 2nd order mid-point rule 
at the interval [r»_i, rj x [0j— 1, 0j] x [</>fc-i, </>&] is written 



1 \- f{rj,0j,<t>K) - f{rj-i,0j,<t>K) rr , x 

4^ ^ An ,l J / / ^/ dxj>S(r,0,<t>) 



J=j-1 K=k-1 



S (r^ _ 1 , ■ 1 , </> fc 1 ) Ar 4 A6jAcf>k- (C5) 



J=i-lK'=fc-l J 
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FIG. 14. Plots for the x-component of the shift /3 X of BBH 
initial data for the case with AH boundary conditions ((36}- 
t]38J> . Parameters in the conditions are chosen as no = 0.1 
(solid red lines), and no = 0.005 (dashed green lines), with 
Q — 0.08 and Q s = 0. Top panel: along the rr-axis. Middle 
panel: along the y-axis. Bottom panel: along the z-axis. 




FIG. 15. Plots for the same model as Fig. [14] but for the 
y-component of the shift f} y . 
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FIG. 16. Plots for the same model as Fig. [14] but for the 
z-component of the shift fi z along the z-axis. 



0.06 
0.04 
0.02 
cd 0.00 
-0.02 
-0.04 
-0.06 





1 




i 


1 


1 








COCP-BH 




























































































l 






— a. = 0.1 - 
n s = o.o - 







y 




4 6 8 



FIG. 17. Plots for the components of the shift # of BBH 
initial data for the case with AH boundary conditions (|36[1 - 
(|38f) . Parameters in the conditions are chosen as no = 0.1 
with a spin parameter Q s — (solid red lines), and Q s = 0.1 
(dashed green lines). The spins are aligned to the orbital 
angular momentum (i.e. parallel to z-axis). Top panel: j3 x 
component along the y-axis. Bottom panel: /3 y component 
along the ic-axis. Solid red curves correspond to those in 
Figs, m \M 
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FIG. 18. Same as Fig |17l but the direction of the spin is 
aligned parallel to the j/-axis. Top panel: /3 X component along 
the z-axis. Middle panel: /3 Z component along the i-axis. 
Bottom panel: fi z component along the y-axis. 




FIG. 19. Surface integral for the exterior and the interior of 
a sphere 



